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ABSTRACT 

The  estimation  of  time  delay  and  Doppler  difference  of  a  signal  arriving  at  two 
physically  separated  sensors  is  investigated  in  this  thesis.  Usually,  modified  cross  power 
spectrum  coupled  with  Doppler  compensation  is  used  to  detect  a  common,  passive  signal 
received  at  two  separated  sensors.  Another  successful  approach  uses  the  cross  coherence 
to  achieve  this  goal.  This  thesis  modifies  these  two  techniques  to  model  the  Doppler 
difference  via  an  autoregressive  (AR)  technique.  Analytical  results  are  derived  and  ex¬ 
perimentally  verified  via  a  computer  simulation.  Performance  at  high  and  low  signal  to 
noise  ratios  (  SNRs  )  is  examined. 
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I.  INTRODUCTION 


This  thesis  investigates  the  use  of  autoregressive  (AR)  models  for  estimating  the 
Doppler  difference  and  differential  time  delay  by  processing  of  a  narrow  band  signal 
emitted  from  a  moving  source  and  received  at  two  physically  separated  sensors.  If  the 
signal  is  received  at  two  different  geographical  positions  even  in  the  presence  of  uncor¬ 
related  noise,  then,  depending  on  the  signal  strength  and  duration,  it  is  possible  to  esti¬ 
mate  the  differential  time  delay. 

Because  the  source  is  moving,  the  signals  that  are  received  at  the  sensors  may  have 
different  frequencies  due  to  the  Doppler  effect.  To  obtain  accurate  differential  time  de¬ 
lay  estimation.  Doppler  difference  compensation  is  usually  required. 

This  compensation  can  be  implemented  by  using  frequency  shifting  of  the  narrow 
band  components  of  the  received  signal.  This  frequency  shift  can  be  obtained  using  a 
Fourier  transform.  In  this  thesis  we  use  an  AR  model  to  detect  the  frequency  shift. 
Using  this  Doppler  compensation,  an  estimate  of  differential  time  delay  can  be  obtained. 
Estimating  the  delay  and  Doppler  using  an  AR  model  can  be  interpreted  as  a  form  of  a 
narrow  band  coherence  procedure,  provided  the  estimations  are  properly  normalized. 

This  thesis  is  arranged  in  six  chapters  and  six  appendices.  Because  the  estimation 
of  the  time  delay  and  Doppler  is  intimately  related  to  the  coherence  between  two  trans¬ 
formed  complex  signals,  an  extensive  investigation  of  coherence  is  given  in  Chapter  II. 
In  Chapter  III,  AR  models,  advantages  of  AR  modeling,  and  AR  model  order  selection 
are  presented.  Chapter  IV  is  devoted  to  the  analysis  and  the  processing  of  noisy  signals 
to  estimate  the  differential  time  delay  and  the  Doppler  difference.  To  estimate  the  dif¬ 
ferential  time  delay,  two  approaches  are  pursued.  AR  model  performance.  Doppler  es¬ 
timation  and  two  types  of  time  delay  estimation  are  examined  in  Chapter  V.  In  the  last 
chapter  some  general  conclusions  of  the  work  carried  out  in  this  thesis  are  presented,  and 
some  suggestions  for  future  investigations  are  given.  Computer  simulation  programs  are 
included  in  Appendices  E  and  F. 
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II.  COHERENCE 


A.  DEFINITION 

The  coefficient  of  coherence  between  two  wide  sense  stationary  random  processes 
is  the  normalized  cross  power  spectral  density  function  defined  by  Wiener  [Ref.  1:  p.  12] 
as 


Vxyif)  = 


GXy(f) 

\  Gxx(f)Gyy(J) 


(2.1) 


where/ denotes  the  frequercy  {Hz)  , 

GX)(f)  is  the  cross  power  spectrum  between  .*(/)  and  y(/)  , 

Gx,{f)  denotes  the  auto  power  spectrum  of  .t(t)  ,  and 
G}i{J)  denotes  the  auto  power  spectrum  of  r(t)  . 

Despite  some  confusion  in  the  literature.  Wiener  intended  for  the  coefficient  of  co¬ 
herence  to  be  complex.  This  is  apparent  since  he  discusses  both  the  modulus  and  the 
argument  of  the  coefficient  of  coherence. 


B.  PROPERTY  OF  THE  COHERENCE  FUNCTION 

The  power  spectral  density  matrix  0(f)  is  positive  semi  definite.  Therefore,  for  two 
random  processes  .v  and  y,  we  see  that 


det  [0(/l]  =  del 


Gxxif )  Gxy(J) 

GyX(f)  Gyy(f) 


(2.2) 


For  real  processes  we  have  </,(/)  =  G'xy{f)  and  thus 

Gxx(J)Gyy(/)-\Gxy(f)\2>0  (2.3) 

where  *  denotes  the  conjugate  of  a  complex  number.  And 

GxxifGyyU)  >  |GJ2  (2.4) 

Furthermore.  G,x{f)  and  Gty(f)  are  nonnegative,  real  functions  of  frequency.  When 
G„{/)  and  GJJ)  are  strictly  positive  definite  (  i.e.,  G,x{f)Gty(f)  >  0  ).  Eq.  (2.4)  can  be  di¬ 
vided  by  GJJ)Gtt(J)  without  changing  the  sense  of  the  inequality. 


This  provides  as  an  upper  bound 


!  yxyi/)  I  v/  (2.5) 

and  since  the  magnitude  of  any  complex  number  is  greater  than  or  equal  to  zero,  we 
have  the  lower  bound 


0  <  I  rxy(f)  I  <  1  (2.6) 

The  magnitude  of  the  coherence  function  is  always  between  zero  and  one.  It  is  zero 
if  the  processes  .v(/)  and  y(i)  are  uncorrelated  and  it  is  equal  to  unity  if  there  exists  a 
linear  relationship  between  x(r)  andy(/)  .  In  order  to  define  the  coherence  it  is  necessary 
that  the  numerator  and  denominator  of  Eq.  (2.1)  are  not  simultaneously  equal  to  zero. 
Coherence  is  not  defined  if  either  auto  spectra  is  zero. 

C.  COHERENCE  ESTIMATION 

If  .Y,(,0  denotes  the  Fast  Fourier  Transform  (FFT)  of  the  /  th  segment  of  jt(m)  at 
frequency  fk  .  then  the  spectral  density  estimates  are  given  by  [Ref.  2:  p.  22] 

v 

4l4)  =  *X|A'^12  (17) 

i=\ 

V 

Gxyifk)  =  (2.8) 

/=! 

V 

4(4)=»Zly'^)>2  (2-9) 

/=! 

where  a  =  • 

.V  =  number  of  segments, 

T  =  segment  length,  and 
ft  =  sampling  frequency. 

Finallv  the  coherence  estimation  is 
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Gxy(fk) 

X  Gxx(fk)Gyy(fk) 


V  Mrtm 


2 


(2.10) 


D.  COHERENCE  OF  NARROW  BAND  SIGNALS  WITH  DIFFERENTIAL  TIME 
DELAY  AND  DIFFERENTIAL  DOPPLER 

The  output  of  the  band  pass  filters  (BPF1  and  BPF2)  in  Figure  1  are  denoted  by 
A' (/l)  and  YUfi)  respectively.  Each  term  represents  the  Fourier  transform  of  the  corre¬ 
sponding  time  series  evaluated  ai  frequency/,  and  time  /.  For  narrow  band  signals  a 
Doppler  shift  corresponds  to  a  frequency  shift.  If  a  signal  arrives  at  the  two  sensors 
having  a  differential  Doppler  shift  as  well  as  a  differential  time  delay,  then  we  see  that 
a  frequency  compensation  and  time  delay  compensation  by  the  appropriate  values  tend 
to  line  up  the  signals  in  frequency  and  time.  This  is  accomplished  by  using  an  additional 
Fourier  transform  in  channel- 1  of  the  processor  and  a  time  delay  in  channel-2  of  the 
processor.  Mathematically,  this  can  be  expressed  as 


C-11) 


Comparing  this  with  Eq.  (2.10),  we  see  that  we  generalized  the  coherence  concept.  We 
also  note  that  the  implementation  resembles  a  correlator,  where  one  of  the  signals  is 
frequency  compensated  and  the  time  delay  corresponds  to  the  delay  operation  of  the 
correlator. 
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Figure  1.  Coherence  estimation  block  diagram. 


If  the  Figure  1  is  redrawn  as  in  Figure  2,  then  it  can  be  interpreted  as  an  FIT 
implementation  as  shown  in  Figure  3. 


Figure  2.  Coherence  estimation  block  diagram  (reinterpretation). 


BPF1 


Figure  3.  Coherence  estimation  block  diagram  using  the  FFT. 


But  this  FFT  has  a  poor  resolution.  To  obtain  a  better  resolution,  an  AR  model 
is  desirable.  This  approach  is  shown  in  Figure  4. 


Figure  4.  Coherence  estimation  block  diagram  using  an  AR  model. 

Figure  1  and  Figure  4  represent  two  difTerent  schemes  to  compute  the  coherence 
function.  Differential  time  delay  and  Doppler  difference  can  be  estimated  by  using  AR 
models  in  the  modeling  of  the  coherence  function. 
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III.  AUTOREGRESSIVE  (AR)  MODELS 
A.  AR  MODELING 

If  the  following  difference  equation  is  satisfied,  the  resulting  structure  is  called  an 
AR  model  of  order  p.  [Ref.  3:  p.  177J 

p 

x(")  =  -  ~  *)  +  u(")  (3.1) 

*=  1 


where  x(n)  —  the  signal  at  instant  n, 
u(/t)  =  the  white  noise  driver,  and 
ak  =  the  k  til  coefficient  of  the  AR  model. 

A  realization  of  the  AR  model  is  illustrated  in  Figure  5. 


Figure  5.  AR  Filter  of  order  p. 
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B.  ADVANTAGES  OF  THE  AR  MODEL 

The  motivation  for  parametric  modeling  of  random  processes  is  the  ability  to  obtain 
better  spectral  estimates  based  upon  the  model  than  estimates  produced  by  classical 
spectral  estimation.  Both  the  periodogram  and  correlation  methods  can  be  used  to  yield 
Power  Spectral  Density  (PSD)  estimates.  The  unavailable  data  or  autocorrelation  se¬ 
quence  (ACS)  values  outside  a  given  window  are  assumed  to  be  zero.  This  kind  of  an 
unrealistic  assumption  leads  to  distortions  in  the  spectral  estimate. 

The  advantages  of  the  AR  approach  are 

1.  AR  spectra  tend  to  have  sharp  peaks,  a  desirable  feature  of  high  resolution  spectral 
estimators. 

2.  AR  parameter  can  be  obtained  as  solutions  to  linear  equations. 

C.  POWER  SPECTRAL  DENSITY 

Eq.  (3.1)  can  be  rewritten  as  follows 

r 

-<•(«) = -  “ A') +  w(") 
k=l 

=  Thku(„-k)  (3.2) 

k=0 

where  h t  =  the  causal  filter  impulse  response. 

Let  us  take  the  Z  -transform  of  Eq(3.2). 

X(r)  =  H(z)U(z)  (3.3) 


Rewriting  Eq.  (3.1)  as 


p 

~k)  =  w(«) 

k= 0 


(3.4) 


where  a0  =  1,  and  taking  the  Z  -transform  gives 

A{z)X{z)  =  l\z) 

Eq.  (3.3)  and  Eq.  (3.5)  can  be  solved  to  obtain  H(z) 


(3.5) 
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(3-6) 


11(2)  = 


1 

A(z) 


and  hence 


X(z)  = 


t'(z) 

A(z) 


(3.7) 


The  Z  -  transform  of  the  output  sequence  (.r(«)}  is  related  to  the  Z  *  transformation 
of  the  input  random  process  u(n)  by  {Ref.  4:  p.  56] 

Pxx(z)  =  X(z)X(  4- )  = - -  PM - “ — j —  (3.8) 

A(z)A  (-V)  A(z)A  (~V) 


The  AR  power  spectral  density  is  obtained  by  substituting  r  =  er-"rr  into  Eq.  (3.8) 
and  scaling  it  by  the  interval  T. 


Par(J)=TP> 


M00l‘ 


Tpa 


e"(J)aaHeF(J) 


(3.9) 


where  c.  =  [  1  .<r!2''T, . e-jwr^H 

=  variance  of  driving  sequence. 

D.  BURG'S  ALGORITHM 

In  practice,  the  autocorrelation  is  usually  not  available,  so  one  must  make  an  AR 
spectral  estimation  based  on  the  available  data  samples.  The  simplest  procedure  to  ob¬ 
tain  an  AR  spectral  estimate  from  data  samples  would  be  to  produce  estimates  of  the 
autocorrelation  sequence  from  the  data.  These  autocorrelation  estimates  would  be  used 
in  lieu  of  the  true  autocorrelation  sequence  in  the  YULE  -  WALKER  equations  to  yield 
the  AR  coefficients.  However  better  results  are  obtained,  particular  for  short  data 
segments,  by  algorithms  that  obtain  the  AR  model  parameters  directly  from  the  data, 
without  explicitly  estimating  the  autocorrelation  function. 

The  Levinson  recursive  solution  to  the  YULE  -  WALKER  equations  relates  the 
pth  order  AR  parameters  to  the  p  -  1  th  order  parameters  as  given  by  [Ref.  3:  p.  211  ] 

ap(")  =  Vi(")  +  kPal-\(P  ~  ")  (3- 10) 
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For  n=  1  to  n  =  p—  l,  the  reflection  coefficients  kF  can  be  found  by  using  the 
known  autocorrelation  function  for  lag  0  to  p  —  1.  So  we  have 


"Xvi  (n)rxxiP~n) 

kp  =  ap(p)  =  " -ppL~ -  (3.11) 

The  recursion  for  the  driving  white  noise  variance  is  given  by 

Pp  =  PjP_1(l  -  lAp|2)  (3.12) 


where  p0  =  r„(0)  . 

But  the  ACS  is  not  available,  hence  we  can  not  calculate  the  reflection  coefficients. 
The  Burg  algorithm  provides  an  estimate  of  the  reflection  coefficients  which  in  turm  are 
obtained  through  a  least  squares  criterion. 

The  forward  linear  prediction  error  is  given  by 


e^n)  =  x(n)  + 


p 

S^j/p{m)x{n  -  m) 

m  =  1 


(3-13) 


while  the  backward  linear  prediction  error  is  given  by 


«£(»)  =  *(«  -  p)  + 


p 

^  dp  ( m)x(n  +  m- p) 


(3.14) 


Substitution  of  Eq.  (3.10)  into  Eq.  (3.13)  and  Eq.  (3.14)  yields  the  recursive  relationships 


dp{n)  =  e£_,(n)  +  V!-i("  ~  1) 

(3.15) 

*£(«)  =  Vi(«  -  1)  +  kp/p^ri) 

(3.16) 

At  each  order  p,  the  arithmetic  mean  of  the  forward  and  backward  linear  prediction 
error  power  (  sample  prediction  error  variance)  is  given  by 
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<■?-+ 


"v"  Z  l*p(w)l2  +  -y  Z  l#«)l 


«=p+i 


(3.17) 


This  expression  is  minimized,  subject  to  the  recursion  given  by  Eq.  (3.15)  and  Eq. 
(3.16).  Thus,  f>f  is  a  function  of  single  parameter,  namely  the  complex  valued  reflection 
coefficient  k„.  Setting  the  complex  derivative  of  Eq.  (3.17)  to  zero 


arf  .  erf 

oRe{kp)  '  1  clm(kp) 


and  solving  for  k .  yields 

w  r  • 


(3.1S) 


v 

-2  ^  4-iMeP-i(n  ~ 

4  =  — - — - - -  (3.19) 

Z  i^-i(w)i2+  Z  i^-i(«-i)2 

«=p+l  /!=/>+ 1 


The  estimation  of  the  reflection  coefficient  represents  the  HARMONIC  mean  be¬ 
tween  the  forward  and  backward  partial  correlation  coefficient,  where 
<?o(w)  =  e$(/i)  =  .r(fl)  and  .V  =  number  of  data  points. 

E.  FINAL  PREDICTION  ERROR  (FPE)  CRITERION 

Because  the  best  choice  of  filter  order  is  not  generally  known  a  priori,  it  is  usually 
necessary  in  practice  to  postulate  several  model  orders.  FPE  is  a  kind  of  criterion  '.ihich 
was  provided  by  Akaike  (Ref.  3:  p.  230].  This  criterion  selects  the  order  of  the  AR 
process  so  that  the  average  error  variance  equals  the  sum  of  the  power  in  the  unpre¬ 
dictable  part  of  the  process  and  of  a  quantity  representing  the  inaccuracies  in  estimating 
the  AR  parameters.  The  FPE  for  an  AR  process  is  defined  as 


FPE{p )  = 


a  N  +  (p  +  1 ) 

Pp  x-{p  +  i) 


(3.20) 
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where  A'  is  the  number  of  data  samples, 
p  is  order  of  the  filter,  and 

Pp  is  the  estimated  white  noise  variance  when  using  a  p  th  order  filter. 
The  order  p  selected  is  the  one  for  which  the  FPE  is  minimum. 
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IV.  DOPPLER  AND  DIFFERENTIAL  TIME  DELAY  ESTIMATION 


A.  DOPPLER  ESTIMATION 

Let  jc(/)  and  >•(/)  be  the  signals  received  by  two  sensors 


x{[)  =  cos{2?r(/' 4-  «j)r} 

(4.1) 

>•(/)  =  cos{2rr  (f  +  a2)t} 

(4.2) 

where/is  the  carrier  frequency, 

a,  ,  arj  are  Doppler  shifts,  and 
/  is  the  time  variable. 

Let  f  be  the  sampling  frequency  which  satisfies  the  Nyquist  theorem.  For  con¬ 
venience.  in  all  derivations  and  simulations  a  sampling  frequency  of  64  Hz  is  used,  to¬ 
gether  with  band  pass  filter  width  of  1  Hz.  We  assume  that  |a,|  <0.5  (  /=  1,2  )  . 
which  implies  that  the  signals  stay  in  the  band  pass  filter  regions  of  their  respective  band 
pass  filters  regardless  of  any  Doppler  shift.  Note,  these  values  can  be  modified  to  arbi¬ 
trary  sampling  rates  and  pass  band  regions.  The  sampling  rate  is  given  by 


fs>2f+\>2{f+  y.j)  (  /=  1,2  )  (4.3) 

The  phase  at  instant  k  and  sampling  time  interval  T  are  given  by 


* 

II 

K» 

^  + 

5? 

(4.4) 

,  2  n{f  +  a2)k 

yVk  ~  j' 

(4.5) 

r=  — 

fs 

(4.6) 

Let  x(«)  denote  the  sampled  analog  signal  x(t)  ( ^„T  then 

x(n  +  m)  =  cos{2  n[f  +  a,) -j-  +  x4>m}  (4.7) 
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(4.8) 


y{n  +  m)  -  cos{2n{f  +  ct2)  y-  +  y<f>m) 

Js 

The  derivations  of  Eq.  (4.7)  and  Eq.  (4.8)  are  given  in  Appendix  A. 

These  signals,  (jc(/i)}  and  {>'(«)},  are  processed  to  detect  the  Doppler  difference  . 
Let  BPF1  be  the  band  pass  filter  centered  at  yj  and  BPF2  be  the  band  pass  filter  centered 
at  f  (where  f  might  equal /I  ).  N  data  points  are  processed  in  the  band  pass  filters,  to 
generate  one  output.  The  inputs  to  BPF1  and  BPF2  at  time  m  are  vectors  XJn)  and 
}’„(«) ,  respectively.  The  size  of  the  input  vector  to  the  filter  is  the  number  of  data  points 
taken  during  1  second  (i.e.,  the  number  of  points  processed  in  the  filter  =  X  —f  ).  The 


input  vectors  are  denoted  by 

Xm(n)  =  [.x(m).  x(m  +  1 ) . .  x{m  +  A'  -  1 )]  (4.9) 

Ym(n)  =  4-  1), . ,y{rn  4-  .V  -  1)]  (4. 10) 


The  filtering  is  performed  using  FFTs  ,  where  successive  FFT  outputs  are  generated  at 
the  input  data  rate.  The  BPF2  output  is  conjugated. 

To  avoid  the  complexity  the  following  four  complex  constant  variables  are  defined. 


.v-i 


Ax  =  Vcos(2*(r  4-  a,) y  }e  ’j2'f  s 

Js 


(4.11) 


,v-i 


Bx  =  X5in^2^ +  /■  ^  S 


n= 0 


(4.12) 


v-i 


Ay  —  ^cos(2 n{f  +  a2)  y-  }e>2^  -v 

n=0  ■'* 


(4.13) 


v-i 


By  =  V  sin{2 n{f  +  a2)  y  }el2~f  s 
Js 


(4.14) 


At  instant  m  the  complex  output  signal  of  BPF1  can  be  calculated  as  follows 
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xm{f)  =  ^cos{2t z{f  +  «,) -jr  +  x4>m)e  Jl'f  v 

«=o  •'* 

.v-i 

=  cos{2tt(/  +  a,)  jr  }  cos  x<f>m  -  sin{2 n(f  +a])~}  sin  x<t>m']e~j2'f— 

n= 0  1  5 

V-I  v-i 

=  cosx4>mY]cos{2r,(f+a])-jr}e'j2z~  -smx4>„,y]s\n{2r,(f+al)~r}€'j2'~\ 
rTn  Js  tZn  Js 


=  ^  COS  -  Bx  sin  x<t>m 

^x$m  _j_  (fix®™  x& n. 

—  ^4  _  Z?r  _  — 


(4.15) 


=  A  V+-iB*  +  A*  ijgi  -j,Q- 


Ax+jBx  aoV±2lLm  ,  Ax~jBx  _MV±1 ±m 
- - - e>  *  fs  + - - - e  j. 


The  complex  signal  XJJ)  contains  two  frequencies  with  different  amplitudes. 

To  understand  the  character  of  A'„,( J)  ,  it  is  important  to  evaluate  which  is  the 
dominant  term  in  the  above  expression. 


V-l 

Ax  +jBx  =  7;  C0S{2 n(f+  «,)  J  }  +j  sin{2 rt(f+  «,)  jr  J]^' v 


n=0 

V-1 


fs 


_  ai)  /*  e  Jl'f  v 

;:=0 


(4.16) 


/7=0 

1  - 
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Ax~jBx  =  ^\cos{2n(f+a:)yr}  -./sin{2jr(/,+ a,)  jjr }]«-  j2r/  x 
*=  0 

-  V 


nssQ 

.V-l 


(4.17j 


-Z-- 


-y:*(2/+  a,)-^- 


/7=0 


I  —  e~ 


\-e-jWf+*i)T 


From  Appendix  C 

I  Ax  +jBx  f  •>  I  -'lx  ~jBx  I  (4. 1 S) 

A,  4-/7?, 

Therefore  — ’—zf — -e>*«  is  the  dominant  term  of  .Vm(/). 

In  similiar  way,  the  complex  output  signal  of  BPF2  can  be  calculated  as  follows 


V-l 


1 «(/)  =  yjCos{2rr(/ 4-  a:)  —  +  v 


>:=0 

\-1 


=  cos{2rlf  4-  x;)  j- }  cosy4>m  -  sin{2n-(/'  +ix2)j-)  sin  x 

>;= 0  *  5 

V-l  V-l 

=  cos  v(j)my^cos{2r.{f  +  <x2)~}e/2'T  -  sin^Y sin{2r(/r4-  x2) y-  }c_y'' ~ 


=  /l,  cos  y4>m  -  By  sin  y<f>„ 


(4.19) 


=  /F 


-5, 


>  *>/ 


-4V  +/5V  •  -/'5V  .  ^ 

=  — •  ■  ■■  ■  -  4.  — ; - —  ^ 

2  2 

4+A  J26^pLm  ,  4-^v 

—  Z  J  f  ~T~  €  Jt 


The  complex  signal  }*(/)  contains  two  frequencies  with  different  amplitudes. 
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Even  though  Y'Jf)  is  similiar  in  form  to  XJJ),  it  is  not  obvious  which  term  is  dom¬ 
inant. 


Ay  + jBy  =  ^  [  cos{2tt(/'+  a2)  ~r  }  +j  sm{2n(f  +  a2)  jr  \\^'f  .v 

/3=0  5 

•V-1 

.  y 


/7=i.) 

v-1 

z 

n=() 

.V-l 

/iasO 

1  _  y2) 

,_r;2.(2/+,2)J- 
1  -  42r*; 


(4.20) 


4  -/A  =  y  [  c°s{2-(/+  «2) -f  }  -j  sin{2 r(f+  «,)  -f  }]42!r/_T 

ts'  A 

,v-i 


/7=U 

.V-l 

-z- 


•  ->_  /7 


1  -  c"72^ 

\-e~J2^T 


(-1-21) 


From  Appendix  C 


l  Ay  ~jBy  |  >  \Ay+jBy\ 


Therefore 


4  ~jB> 


■  eycm 


is  the  dominant  term  of  }^(/). 


(4.22) 
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Since  the  two  output  sequences  [XJJ)}  and  { ¥'„{/)}  have  two  sinusoidal  components 
each,  their  product  {A ’„(/)!"(/)}  has  four  sinusoidal  components.  The  four  frequencies 
are 

„  0*i -<*2)  ,  (*2-«i)  „  (2/+a,+oc2)  „  (2/  +  a,  4-  a2) 

a)t  =  2tt - - -  ,co2  =  2n - - -  ,  cu3  =  27r - - -  ,  cod  =  — 2tt - - - 

Js  Js  Js  Js 

with  co,  the  sinusoid  with  the  largest  amplitude  and  co2  the  sinusoid  with  the  smallest 
amplitude. 

The  product  of  the  output  of  the  filters  is  given  by 


u/>  rim  =  ^4^  A-A  £fz 


(Ax  +jBr)(Av  -jBv) 


2  ■b^"(aii  a2’  x^m- yfim) 


where  Fla,.  ,<+>„,  y<j>J  represents  the  three  low  amplitude  frequency  terms. 

When  using  the  AR  model  as  described  in  chapter  III.  at  a  high  SXR  (  i.e., 
SXR  ->  00  )  ,  a  4th  order  AR  model  detects  all  of  the  four  frequencies  given  above. 
When  the  SXR  is  low  (i.e.,  20  dB  ),  only  one  dominant  sinusoid  is  detected  with  fre¬ 
quency 


Since  we  can  detect  co,  and  know  the  value  of/  ,  a,  -  a2  can  be  estimated  by  using  an 
AR  model. 


B.  DIFFERENTIAL  TIME  DELAY  ESTIMATION 

To  detect  and  localize  a  target,  it  is  important  to  find  the  differential  time  delay  of 
signals  arriving  at  two  sensors.  Let  us  assume  that  signals  at  the  two  sensors  are  as 
given  in  Figure  6  (  i.e..  zero  differential  time  delay  )  For  the  remainder  of  the  figures,  the 
time  avis  is  scaled  to  be  64  points  per  second. 
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Figure  6.  Receiving  signals  at  sensor  1.2  (SNR=  100,  no  differential  delay). 


We  can  estimate  the  differential  time  delay  and  differential  Doppler  using  two  dif¬ 
ferent  approaches.  The  first  approach  uses  the  cross  power  spectrum  while  the  second 
approach  uses  the  coherence. 

1.  Differential  time  delay  and  differential  Doppler  estimation  using  the  cross  poner 
spectrum. 

Figure  4  shows  the  block  diagram  of  the  coherence  estimator  using  AR  mod¬ 
eling.  In  this  figure,  the  output  of  the  FFT  can  be  interpreted  as  the  cross  power  spec¬ 
trum.  Using  this  cross  power  spectrum,  we  can  estimate  the  differential  time  delay  and 
the  differential  Doppler.  But  when  we  use  the  highest  peak  of  the  cross  power  spectrum, 
the  peak  is  somtimes  not  detected  at  the  proper  time  delay  nor  at  the  proper  Doppler 
shift.  We  will  show  a  special  case  in  which  the  peak  of  the  power  spectrum  is  not  de¬ 
tected  at  the  proper  time  delay  nor  at  the  proper  Doppler  values. 
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a.  Special  case. 

For  our  test  case  the  data  duration  is  six  seconds,  two  seconds  of  the  noisy 
signal  and  four  seconds  of  noise  only.  Linear  transformation  of  this  data  leads  to  one 
of  three  types  of  outputs.  The  first  type  represents  full  information,  while  the  second 
type  represents  partial  information.  The  third  type  represents  a  noise  only  condition. 
Generally  when  two  signals  are  lined  up  in  time,  the  AR  model  should  give  the  highest 
output  power.  But  this  is  wrong  in  some  cases.  In  the  full  information  case,  the  mag¬ 
nitude  of  transformed  signal  is  high  and  constant.  For  some  reduced  information  case, 
even  though  the  magnitude  of  the  transformed  signal  decreases,  it  still  may  have  large 
magnitude  of  spectral  components.  This  phenomenon  can  be  explained  as  follows. 

Define  two  functions  depending  on  k 


X—l—k 


t Ay  =  Yj  c°s(2^/2  -jr  )ei2'f  A 


n=0 


A'-l-fc 

kBv=  y  sin(27r/2  —  )el2~f~ 
n- 0  Js 


(4.25) 

(4.26) 


where  f2  =/+  a2. 

A  denotes  the  number  of  lost  data  at  BPF2  input  vector. 


We  assume  that  we  know  when  the  BPF1  signal  starts.  If  this  information 
is  not  available,  we  need  to  examine  the  signal  at  the  output  of  BPF1.  to  obtain  a  can¬ 
didate  time  frame.  Consider  the  input  of  AR  model  as  shown  in  Figure  7. 

The  input  data  size  to  the  AR  model  is  the  number  of  linearly  transformed 
data  points  during  a  given  period  (i.e.,  input  size  is  A'  =f  ).  The  BPF1  output  represents 
full  information  (i.e.,  each  output  element  is  produced  from  the  signal  information 
points).  Therefore,  as  shown  in  the  previous  section,  the  dominant  output  sequence  is 

A/0,  A,<y\  A/>2,  kxe>’*\ .  A/^A'-'  (4.27) 


where 


A,  = 


1  1  - 

2  n  1 


(4.28) 
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Figure  7.  AR  model  and  its  driver  source. 

(I )  0  second  delay  case.  We  already  defined  the  input  vector  Xm(n)  and 

}'„(//)  to  BPF1  and  BPF2.  This  vector  can  be  interpreted  as  a  time  snap  shot.  Fach 
output  of  the  BPF  requires  A'  input  points.  If  we  require  A'  output  points  from  the  BPF, 
even  using  maximum  overlap  in  the  processing,  we  require  at  least  2Ar  -  1  points  at  the 
BPF  input.  Figure  8  shows  the  2N  input  to  BPF1  and  BPF2.  Both  2 A'  input  points 
contain  the  signal  and  some  noise.  As  shown  in  Figure  8,  in  the  zero  delay  case,  the  two 
input  vectors  Xm{n)  and  Ym{n)  (  0  <,  m  <  A'  )  provide  full  information.  So  the  lin¬ 
early  transformed  output  X„{ft)  and  r„(fk)  also  represent  full  information. 

The  dominant  part  of  the  BPF2  output  is  the  sequence 

*2e~yo,  k2e~J>*\  k2e~J^2,  k2e~J^3 . .  k 2e~J><t‘lf-'  (4.29) 

where 


k2~  2 


1  l-e 


1  —  e 


(4.30) 


21 


Figure  8.  Two  signals,  0  second  delay  at  the  sensors. 

The  dominant  input  sequence  to  the  AR  model  is  given  by 


,  A  =  0,1,2 . .V—  1  }  (4.31) 


where  the  magnitude  is  klk2 

f2j  -1  second  delay  case.  Figure  9  shows  the  2.Y  input  to  BPF1  and 
BPF2.  The  2.Y  input  data  points  of  BPF1  contain  the  signal  and  noi^e.  During  the  first 
N  points  the  input  to  BPF2  contains  the  signal  plus  noise  while  during  the  second  .Y 
points  only  noise  is  present.  The  input  vector  of  BPF1  contains  full  information, 
but  input  vector  Fm(«)  of  BPF2  does  not  contain  full  information.  In  the  }'„(«)  case, 
every  element  of  F0(/7)  is  full  information.  But  when  k  is  not  equal  to  zero  then  )\(n) 
consists  of  A*  —  k  information  elements  and  k  noise  elements. 
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The  output  from  BPF2  represents  partial  information  while  the  out¬ 
put  from  BPF1  represents  full  information.  The  BPF2  output  sequence  is  given  by 


{  kA±jJl"B±.  g/A  +  Jdl  JhEl L  g-/A  ^  _  0,1, .--A  —  1  } 


(4.32) 


The  derivation  of  Eq.  (4.32)  is  given  in  Appendix  B. 


kAyt—jBy 

Assuming  that - ^ is  the  dominant  term,  which  is  a  valid  assumption  for  our 

test  case  with  f  =  64/*,  /  =  23/*,  a,  =  0.4/*  and  a2  =  0.001  Hz,  then  the  input  se¬ 
quence  to  the  AR  model  is  dominated  by 


{* 


,  <**>-*)'  k  =  0X2 . A'  — 


1} 


(4.33) 


with  ki  is  defined  earlier  (  Eq.  (4.2S)  ). 
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The  magnitude  of  dominant  input  sequence  is  obtained  by  examining 


<v-i  -k 


kAy  -jkBy  =  V  [  cos{2 n{f  +  a2)}  -j  %\n{2n{f  +  a2)}]e/2'/  ,v 


n=0 

X-l-k 


■l 


n= 0 
:V-1  -k 


(4.34) 


■Z 


g-J2"*2  v 


n=0 

n  x~k 

1  _ 


When  we  compare  the  magnitude  of  the  sequence  in  Eq.  (4.27)  and  the  magnitude  of  the 
dominant  term  in  Eq.  (4.32),  the  magnitude  of  \k2\  >  \kAi.—jiBi\l2  where 


1  1  -  e~J2rx'- 

2  2  i 

1  —  e~‘  2  ~ 


A-.  = 


and 

■  „  n-  s~k 

k-{v  -Jk»Y  1  1  —  e  7  v  a 2 

-)  =  ~T  1 

\~e  J  .v  (4.35) 

_  »  **  1  ^ 

2  1  -e-^jr  2  1 -e“7'2'aj"v 

We  note  that  Eq.  (4.35)  has  two  frequencies  (  i.e.,  0  and  —  2n  ). 

The  magnitudes  of  both  terms  are - - - — 

\l-e-;2n,2jr\ 

Now,  note  that  as  long  as  |a2|  <  -jr ,  then  the  magnitudes  of  both  terms  in  Eq.  (4.35) 
are  larger  than  I  k2  \ 


Due  to  the  symmetry  a  positive  delay  results  in  a  similar  reponse  to  a  negative  delay,  and 
hence  Eq.  (4.35)  holds  also  for  delay  of  +  l  second.  Equation  (4.35)  can  be  interpreted 
as  having  large  responses  at  shifted  Doppler  values  at  a  delay  of  1  second  and  a  delay 


of — 1  second.  Proper  delay  can  not  be  established  nor  can  the  Doppler  value  be  estab¬ 
lished.  This  cross  power  spectrum  algorithm  will  not  give  good  information  about  either 
Doppler  difference  or  time  delay. 

b.  Modified  cross  power  spectrum. 

To  detect  the  signal,  input  signals  to  the  AR  model  should  be  preprocessed 
as  shown  Figure  10.  If  the  BPF2  output  signal  is  magnitude  normalized,  then  this  nor¬ 
malized  signal  retains  good  phase  information.  When  we  normalize  the  BPF2  output, 
one  of  tu:o  conditions  can  occur.  The  first  condition  (  i.e.,  synchronized  )  leads  to  the 
magnitude  normalization  of  Eq.  (4.29)  and  provides  accurate  frequency  and  delay  esti¬ 
mation.  The  second  condition  occurrs  when  the  signal  in  channel-2  is  not  lined  up  in 
time  with  the  signal  in  channel-1  (  i.e.,  during  the  delay  search  ).  Under  this  condition 
smaller  peaks  in  the  time  Doppler  plane  appear  at  incorrect  values  of  time  delay  and 
Doppler.  Above  70  dB  SNR,  the  correct  peak  is  the  dominant  one  and  allows  proper 
estimation.  Information  from  contour  plots  can  be  used  at  values  of  SNR  between  70 
and  10  dB  (  see  Appendix  D  ). 


Figure  10.  Modified  cross  power  spectrum  block  diagram. 


2.  Differential  time  delay  and  differential  Doppler  estimation  using  the  coherence. 
Another  w’ay  detecting  the  signal  is  normalization  of  the  AR  output  with  the 
squared  root  of  the  product  of  P„{f)  (  auto  spectrum  of  x  )  and  P„(f)  (auto  spectrum  of 
y  ).  This  can  be  done  using  two  additional  AR  models  as  shown  in  Figure  11. 
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Figure  11.  AR  model  of  coherence. 
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Using  the  AR  approach,  ,  p,)Xf  ,  pyt  ,  and  AR  coefficients  can  be  ob¬ 
tained.  Therefore  we  get 


PJJ)  =  - 


TP* 


M«(/)l 2 


TPyy 

Pyy(J)  =  - ~- 

"  \Ayy{/)\2 


Tp 


xyxy 


Pxyxyif)  2 

I  ^xyxyif)  I 


(4.36) 


(4.37) 


(4.38) 
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(4.37) 


Py>V)  = 


TPyy 

\Ayy(f)\2 


Pxyxyif)  = 


Ppxyxy 

MW/)  I2 


(4.38) 


where  P„(/)  is  the  power  spectrum  of  A", , 

Pyy{f)  is  the  power  spectrum  of  T,.d , 

PtyJ f)  is  the  power  spectrum  of  X,  Y]_d  , 

AXX{J)  is  the  AR  power  spectral  density  of  X, , 

Ayy{j)  is  the  AR  power  spectral  density  of  Y]_d  , 

Axyxy(J)  is  the  AR  power  spectral  density  of  X,Yt-d  cross  term. 

p„  is  the  driving  noise  variance  of  A’, , 

pyy  is  the  driving  noise  variance  of  Y,.d  , 

piv„  is  the  drivins  noise  variance  of  X,Y]  d  •  and 

r-+. 

Js 


We  obtain  a  coherence  estimate  from  the  three  power  spectra,  by  assuming  that 
I PJJ)  1 |  P„m(J)  I  . 


IPxyifll2 

\PxxV)Pyy[f)l 

1  P 

xv  xy  m  \  Pxyxy  \Axx(J)Ayy(J)\2 

1  Pxxif)Pyy<J)  I  f  PxxPyy  I  Axyxy{f)  I 


A 

Y  = 


Pxvxv 

PxxPyy 


\Axx{J)AyM‘ 


Axyxyif)  I 


(4.39) 


(4.40) 


This  approach  provides  good  information  about  the  differential  time  delay  and 
differential  Doppler  down  to  SXRs  of  20  dB. 
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V.  RESULTS 


This  chapter  presents  graphical  results  obtained  by  applying  the  analytical  results 
of  the  previous  chapters  to  specific  examples  and  carry  ing  out  the  required  computations 
on  the  computer. 

A.  AR  MODEL 

This  section  provides  the  AR  model  performance  tests.  As  we  discussed  in  chapter 
III,  Burg's  algorithm  is  applied  to  this  AR  model.  The  FPE  criterion  is  use  to  find  the 
AR  model  order.  Figure  12  is  the  power  spectrum  of  two  test  signals  sin(<x>,7«)  and 
cos(wj7Vi)  .  Figure  13  shows  the  power  spectrum  of  two  other  test  signals  sin {(o2Tn)  and 
cos [u)2Tn).  Both  sets  of  test  signals  in  these  figures  have  an  S.XR  of  20  dB.  The  fre¬ 
quency/  is  13.45  II:  and /  is  2 '  .45  Hz  and  T  is  -jtj-  second. 

Theoretically,  the  power  spectrum  of  these  signals  have  two  impulse  functions  at/ 
and  — /  (  /  =  1  ,2  ).  Figures  12  and  13  show  experimental  results  agreeing  with  the 
theoretical  results.  Figures  12  and  13  indicate  which  frequency  has  high  power  but  they 
provide  no  phase  information. 

The  AR  model  is  sensitive  to  frequency  but  is  not  sensitive  to  the  phase.  For  any 
selection  of  frequency/  which  satisfies  the  Nyquist  theorem,  the  AR  model  provides 
good  frequency  estimation  provided  the  SHR  is  sufficient  large. 
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Figure  12.  AR  model  performance  test  1  (SNR  =  20  dB). 
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INPUT  SIN(wTN)  F=23.45  T-1/64 
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Figure  13.  AR  model  performance  test  2  (SNR  =  20  dB). 
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B.  DOPPLER  ESTIMATION 

In  section  A  of  chapter  IV,  we  discussed  five  items. 

1.  BPF1  and  BPF2  outputs  have  two  sinusoidal  components  each. 

2.  For  BPF1.  the  dominant  frequency  is  designated  as  and  for  BPF2.  the  dominant 
frequency  is  designated  as  a>y  where 

2ti  if  +  a,) 

= - - - 


—2k(/  +  a2) 


3.  The  power  spectrum  of  cross  term  of  BPF1  and  BPF2  has  four  components  which 
can  be  detected  in  very  high  SXR  (i.e.  SXR  =  100  dB). 

4.  When  the  S.XR  is  low  (i.e.  SXR  =  20  dB).  only  one  dominant  frequency  compo¬ 
nent  is  detected. 

5.  The  dominant  frequency  of  the  power  spectrum  Pxyxi{f)  corresponds  to  the  differ¬ 
ential  Doppler  frequency. 

Figure  14  shows  the  power  spectrum  of  the  BPF1  output  when  /  is  13  Hz.  a,  is 
0.45  Hz.  and  /  is  64  Hz.  Figure  15  shows  the  power  spectrum  of  the  BPF2  output 
when  /is  13  Hz.  a,  is  —0.45  Hz.  and  f  is  64  Hz.  As  expected  both  figures  show  two 
spectral  lines  when  the  SXR  is  100  dB.  The  dominant  frequencies  are  located  at /equal 
13.45  Hz  for  BPF1.  and  at /equal  —12.55  Hz  for  BPF2.  Figure  16  shows  the  power 
spectrum  of  the  cross  term  of  BPF1  and  BPF2.  This  figure  shows  four  spectral  lines  at 
an  SXR  of  100  dB.  while  only  one  dominant  frequency  is  detected  at  an  S.XR  of  20  dB. 
Dominant  frequencies  are  always  located  between  —1  Hz  and  I  Hz.  Figure  17  shows 
a  subplot  of  Figure  16  for  frequency  between  -1  Hz  and  1  Hz.  This  figure  shows  the 
dominant  frequency/  =  .9  Hz  and  the  smaller  spectral  component  at /^  -.9  Hz  for  an 
SXR  of  100  dB.  Only  one  dominant  frequency  can  be  detected  at / equal  0.9  Hz  for  an 
SXR  of  20  dB.  Doppler  difference  from  the  dominant  spectral  component  corresponds 
to  the  true  Doppler  difference  a,  —  a2  which  is  0.9  Hz 

Figure  IS  shows  a  comparison  for  SXR  of  20  dB  and  10  dB.  Figure  19  is  subplot 
of  Figure  IS  for  — 1  Hz</<  1  Hz.  From  this  figure,  when  SXR  =  10  dB,  the  dominant 
frequency  is  shifted  a  little  bit  from  the  true  Doppler  difference.  We  have  shown  the  five 
issues  as  we  have  discussed  in  the  previous  chapter.  In  summary,  for  SXRs  greater  than 
or  equal  to  20  dB.  the  power  spectrum  of  the  cross  term  provides  good  Doppler  esti¬ 
mation. 
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Figure  14.  Power  spectrum  of  the  BPF1  output. 
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Figure  15.  Power  spectrum  of  the  BPF2  output. 
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Figure  17.  Subplot  of  the  cross  poner  spectrum. 
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Figure  18.  .  Comparison  the  Doppler  estimation  at  two  different  SNR's. 
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Figure  19.  Subplot  of  the  Doppler  estimation. 
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C.  DIFFERENTIAL  TIME  DELAY  ESTIMATION 

In  Chapter  IV  section  B.l.a  we  discussed  the  problems  associated  with  using  un¬ 
normalized  quantities  and  suggested  two  possible  approaches  to  solve  the  problem.  In 
this  section  we  will  show  that  for  the  case  of  |a,|  <--  Hz  ,  the  power  spectrum  of  a 
nondelayed  signal  is  smaller  than  the  power  spectrum  of  a  signal  delayed  by  one  second 
(  i.e.,  see  discussion  following  Eq.  (4.35)  ). 

To  estimate  the  time  delay,  we  suggest  the  two  different  approaches.  The  first  ap¬ 
proach  is  the  normalization  of  the  BPF2  output.  The  second  approach  normalizes  the 
cross  term  ?,„,.(/)  with  the  power  spectrum  of  BPF1  P„(/)  and  the  power  spectrum  of 
BPF2  Pn{f)  through  the  equation  defined  by  Eq.  (4.40). 

The  following  figures  are  simulations  for  /=  23  Hz,  a,  =  0.23  Hz,  a2  =  —0.02  Hz, 

T  —  second,  delav  of  0  and  an  SXR  of  100  dB.  Ficure  20  shows  the  surface  plot  of 

04 

P,yMi(f  while  Figure  21  shows  the  contour  plot  of  Figure  22  shows  the  cross 

section  plot  of  Figure  20  with  respect  to  the  DELAY  axis.  This  figure  shows  that  the 
power  spectrum  has  a  spurious  peak  at  delay  of  one  second.  At  the  proper  delay  (i.e., 
0  second  delay)  P,m(J)  has  a  relatively  small  value. 


P,m(f)  is  affected  by  two  terms,  namely  and 


\A,Uf)  \2 


.  Fieure  23  shows  the 


maximum  power  spectrum  term  of  the  transfer  function  and  Figure  24  shows  the  driving 
noise  variance.  At  a  delay  of  ±  1  second  and  of  0  second  the  transfer  function  shows  a 
peak  while  the  variance  of  driving  noise  has  small  value.  In  this  case,  P ,>„(/),  the  product 
of  the  form  pu  ^  \/)p  at  ^  secon^  delay  has  a  smaller  value  than  that  at  any  other 
delay  location.  Using  the  scheme  shown  in  Figure  10,  we  normalize  the  BPF2  output 
and  present  the  result  in  Figure  25.  We  can  detect  the  proper  time  delay  and  the  proper 
Doppler  difference.  Figures  26  and  27  show  the  corresponding  contour  plot  and  the 
power  spectrum  of  the  normalized  cross  term.  For  any  value  of  time  delay  and  Doppler 
difference,  this  approach  provides  good  information  about  differential  Doppler  and  dif¬ 
ferential  time  delay  provided  the  SXR  is  greater  than  or  equal  to  70  dB. 

Results  of  the  second  approach  (  Figure  11  )  are  demonstrated  in  Figures  28  and 
29.  This  approach  also  gives  good  information  about  differential  time  delay  and  differ¬ 
ential  Doppler.  This  approach  provides  good  time  delay  and  Doppler  estimation  for 
SXR  greater  than  or  equal  to  20  dB. 

If  we  use  the  contour  shape,  we  can  also  estimate  the  differential  time  delay  and 
differential  Doppler  for  SXR  of  less  than  20  dB.  Differential  time  delay  and  differential 
Doppler  estimation  can  be  extended  down  to  SXRs  of  about  0  dB  when  using  the  con¬ 
tour  of  the  coherence  surface.  Figure  30  is  a  contour  plot  of  using  this  second  approach 
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( i.e..  coherence  approach  )  for  at  =  0.23  Hz,  a2  =  -  0.02  Hz,  Delay  =  0,  and  an  SXR  of 
0  dB.  The  contour  shape  of  the  second  approach  has  the  shape  of  the  letter  X.  When 
the  SXR  is  greater  or  equal  to  0  dB,  the  proper  time  delay  and  Doppler  difference  seems 
to  be  located  at  the  cross  over  point  of  the  X.  When  the  SXR  is  less  than  0  dB.  the 
contour  shape  does  not  shows  an  X  clearly.  In  the  coherence  approach,  Figure  31 
provides  the  contour  plot,  with  channel- 1  containing  only  noise  and  channel-2  contain¬ 
ing  signal  plus  noise  (  SXR  =  0  dB  ).  Figure  32  shows  the  contour  plot,  when  channel-2 
contains  only  noise  and  channel- 1  contains  signal  plus  noise  (  SXR  =  0  dB  ). 
Figure  33  shows  the  contour  plot,  when  both  channels  contain  noise  only.  The  above 
three  figures  do  not  show  an  X  clearly,  hence  do  not  allow  signal  detection  for  esti¬ 
mation  of  any  parameter.  Using  this  information,  we  can  distinguish  the  signal  combi¬ 
nations  from  the  noise  combinations. 

When  we  estimate  the  time  delay  and  Doppler  difference  using  the  contour  plot,  the 
modified  cross  power  spectral  approach  has  some  problems  as  discussed  in 
Appendix  D. 
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Figure  20.  Surface  plot  of  the  power  spectrum. 
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Figure  22.  Maximum  power  spectrum. 
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Figure  23.  Maximum  power  spectrum  of  the  transfer  function. 
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Figure 


25.  Maximum  modified  cross  power  spectrum  of  the  transfer  function. 
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Figure  29.  Contour  plot  of  the  coherence. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  thesis,  the  differential  time  delay  and  differential  Doppler  are  estimated  using 
two  different  approaches.  The  first  approach  uses  a  modified  cross  power  spectrum  al¬ 
gorithm  while  the  second  approach  uses  a  coherence  algorithm.  In  both  cases 
autoregressive  modeling  of  the  required  cross  spectral  component  is  used.  The  differen¬ 
tial  time  delay  and  Doppler  difference  can  be  estimated  using  a  threshold  at  high  SXRs 
and  a  contour  plot  at  low  SXRs. 

At  high  S.XRs  (  i.e.,  SXR  >  70  dB  ),  the  first  approach  locates  the  dominant  peak 
at  the  proper  time  delay  and  the  proper  Doppler.  The  second  approach  utilizes  two 
additional  AR  models  and  two  additional  FFTs  to  obtain  the  autopower  spectrum  of 
channel- 1  and  channel-2.  At  moderate  S.XRs  ( i.e.,  SXR  >  20  dB  ).  the  second  approach 
has  the  highest  coherence  peak  at  the  proper  time  delay  and  Doppler  frequency.  When 
the  SXR  >  0  dB.  the  contour  plot  of  the  coherence  has  the  shape  of  the  letter  X.  The 
differential  time  delay  and  Doppler  difference  can  be  located  at  the  cross  over  point  of 
the  X.  In  high  S.XRs  we  can  use  the  highest  peak  to  detect  the  time  delay  and  Doppler 
frequency  using  either  approach.  When  we  use  the  contour  plot,  the  information  about 
time  delay  and  Doppler  depends  on  the  location  and  the  form  of  the  cross  over  point 
of  the  X. 

We  can  estimate  the  differential  time  delay  and  differential  Doppler  using  the  AR 
modeling  of  coherence.  But  we  can  not  get  the  numerically  correct  value  of  the  coher¬ 
ence  coefficient.  The  following  three  points  are  recommended  for  future  study. 

1.  To  get  a  cross  power  spectrum,  we  assumed  that  I  P^if) 1  ~  I  P,JJ)  1 2  ■  This  was 
required  since  we  can  not  get  a  cross  power  spectrum  using  an  AR  model.  Im¬ 
provement  of  the  cross  power  spectrum  estimates  using  an  AR  model  can  lead  to 
improvement  of  the  algorithm. 

2.  Our  first  approach  is  not  normalized  by  the  auto  power  spectra.  If  we  find  the 
corresponding  auto  power  spectra,  we  can  also  obtain  a  coherence  coefficient  esti¬ 
mate. 

3.  To  get  a  better  information  at  a  low  SXR,  we  need  to  study  the  characteristics  of 
the  contour  plots. 
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APPENDIX  A.  PHASE  DERIVATION 


Define  the  phase  at  instant  k. 


X<t>k  = 


y^k  ~ 


Injf+txjk 

fs 

2n [f  4-  a2)k 

fs 


Let  *(«)  denote  the  sampled  analog  signal  -x(0 1  r=„r 
then 

x(n)  =  cos{2n(f  +  a,) 7- } 

Js 

=  cos{2tt(/  +  a,)  j-  +  x<t>0] 


x(n  +  1)  =  cos{2 n(f  +  a;)  n  jr  -  } 

=  cos{2 n[f  +  «j)  y-  +  ■"  ^  } 

/j 

=  cos{27t(/'+a1)-^-  +  JC0,} 


.v(«  +  2)  =  cos{2rr(/*  +  a,)  } 

JS 

=  cos{2  n{f  +  a,)  7-  +  —  } 

Js  Js 

=  cos{2?r(/'  +  +  x<j)2} 

Js 


x(n  4-  m)  =  cos{27t(/'  +  a;)  - 7-^- } 

Js 

w  27r(/'+a,)m 

=  cos{2tt(/  +  a,)  —  + - - - } 

•/  T  /j 

“  cos{2r:(/'  +  a,)  jr  +  >m) 


(4.4) 


(4.5) 


(a.l) 


(fl-2) 


(a-3) 


(4.7) 


Similiarily 
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A. 


APPENDIX  B.  OUTPUT  OF  BPF2 

Define  two  functions  depending  on  k 

X-]-k 

kAy=  Y  cos(27t/2  7"  (*^-25) 

n= 0  As 

A  —  1  —k 

kBy=  Y  sin(2^/2  y-  )e/2c/~  (4.26) 

«-o 

where  =/+  ct2. 

FULL  INFORMATION  CASE  (  N  POINTS  ) 

A'  points  are  data  points  at  BPF2  input. 


Y-l 


>  wCO  =  V  cos(2-/:  ~  +  y4>,y^  V 
>1=0  ' 


Y-l 


=  cos(2,t/2  -J- )  cos y4>m  -  sin(2^/,  y- )  sin  y^m}^2"f  -v 

/:=0  5 


A  — 1 


.Y-l 


=  cos  v0m y^cos(2-/2  ~  )eilz  s  -  sin  V  sin(2)r/2  y  )<?  y2“  \ 


«=0 


/ 


>1=0 


(6.1) 


=  Q.4y  cos  v(/>m  -  sin  y4>m 

£^2  +  g^-*0"7  D  -  e~ly<t'" 

—  0 -4y  2  2 / 

oAv  +j0Bv  0Av-j0Bv  _j6 
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B.  PARTIAL  INFORMATION  CASE  (  N-l  POINTS  ) 

.V  —  1  points  are  data  points  at  BPF2  input. 

.v-1-1 

y'mif)  =  X  cos(2  rtf2  f  +  y4>J^fJk 
>1=0  'fs 
A*— 1  —  1 

=  ^  {  cos(2 re/2  -j- )  cos>.0m  -  sin(27r/2  jr  )  sin  y4>m}e,'2'f~ 

n= 0  ^  5 

A  — 1— 1  A— 1—1 

~  cos  y<t>m  Yj  cos(2r/2  Y  )^2'  ~  -  sin  y4>m  Y  sin(27r/2  y- )e_/2'~  ^<2) 

n—0  n= 0 

=  v4y  cos  y4>m  -  ]By  sin  y<j>m 

.  e>>*~  -f  n  e1**”  -  e~jy6" 

=  ,4  2  - , By  - 

i  Av  +j]  K  .  \Av-j\Bv 

=  — — — -  +  — — — -  e  J>Pn 


PARTIAL  INFORMATION  CASE  (  N-K  POINTS  ) 

A"  —  k  points  are  data  points  at  BPF2  input. 

X—}—k 

>',*,(/)  =  7  cos(2tt/2  ~  +  v<p„y‘'J~ 

„  A 


n= u 


=  ^  {  cos(2r/2  Y  )  cos  y<t>m  -  sin{2nf2  y- )  sin  y<bm}^f  s 

X-\  -k  X—]—k 

=  C0Sv4m  >  cos(2t:/2  — y  '•  V  -  sin y4>m  >  sin(27t/2  —  )tTy‘'T  ^ 

,._n  Js  „_n  •'* 


—  kAy  cos  y<pm  —  kBy  sin  y4>m 

<jy°»  +  fjy**  _  g-Jytn 


k^y  ^ 


- T} 


kAV  +jkBv  .  ^  kAy  j](By  . 

- ; - —  <*  j - ; - — 

2  2 
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APPENDIX  C.  AMPLITUDE  COMPARISON 


Let  C,  and  C2  be  the  complex  amplitude. 

1  -  e 127:3 


Ci- 


c2  = 


I  -  e-nn* 

1  _  e-jl-0-J+*)  -Jf 


(cl) 


(c.2) 


where  |  a  |  <  0.5,  and  N  >  2f+  1  >  2(f  +  a).  The  numerators  of  |  C,  I  and  |  C2 1  are  same, 
because  the  magnitude  of  the  complex  conjugate  is  the  same.  If  we  show  that 
1 1  —  "^1  <  1 1  —  |  ,  then  we  can  say  I  C,  |  >  |  C2| .  As  shown  in  Figure  34 

1 1  -  =  |,.i/|  and  1 1  -  =  I  Dl\  . 


Figure  34.  Magnitudes  of  two  complex  number. 


The  smallest  angle  LlOB  is 
For  positive  / 


2n(2f  +  a)  2n{X—2f'—a) 

- — - or - r- - 


2rr(l  —  a) 
A' 


LlOB  = 


2n(2f  +  a) 

2na. 

X 

> 

X 

LlOA 


(c.3) 


59 


APPENDIX  D.  CONTOUR  PLOTS  OF  MODIFIED  CROSS  POWER 

SPECTRUM 


For  SXRs  greater  than  or  equal  to  10  dB  the  contour  plot  of  the  modified  cross 
power  spectrum  algorithm  has  the  shape  of  a  boomerang.  The  shape  of  the  boomerang 
depends  on  a2.  Figure  35  and  Figure  36  show  the  contour  plot  of  the  modified  cross 
power  spectrum.  The  two  figures  show  the  boomerang  shape  and  their  different  di¬ 
rections.  For  positive  a2  (  case  1  )  and  negative  a2  (  case  2  )  the  opening  of  the 
boomerang  is  toward  the  left  and  right,  respectively.  In  the  both  cases,  we  can  estimate 
the  time  delay  and  Doppler  frequency  using  these  contour  plots.  The  proper  time  delay 
and  Doppler  difference  can  be  detected  by  locating  the  point  of  symmetry  of  the 
boomerang.  The  Figure  37  shows  a  contour  plot  when  a2  is  —0.02  Hz.  Since  a,  is  very 
small,  the  contour  plot  looks  like  a  straight  line.  In  this  case,  the  differential  time  delay 
and  Doppler  difference  can  still  be  estimated  by  locating  the  point  of  symmetry  of  the 
boomerang  (  i.e.,  center  point  ).  Figure  3S  shows  the  contour  plot,  when  channel- 1 
contains  noise  only  and  channel-2  contains  signal  plus  noise  (  i.e.,  SSR  =  10  dB  ).  Fig¬ 
ure  39  shows  the  contour  plot,  when  channel-2  contains  noise  only  and  channel- 1  con¬ 
tains  signal  and  noise  (  i.e.,  SSR  =  10  dB  ).  Figure  40  shows  the  contour  plot,  when 
both  channels  contain  noise  only.  Using  the  modified  cross  power  spectrum  technique, 
the  three  types  of  noise  only  set  ups  appear  as  sinnliar  plots  (  see  Figures  3S.  39  and 
40  ).  Figure  3S  seems  to  indicate  the  presence  of  signals  in  both  channels,  while  Figures 
39  and  40  appear  more  like  noise.  Hence,  contour  plots  of  the  modified  cross  power 
spectrum  do  not  allow  clear  identification.  In  this  case,  it  will  be  difficult  to  distinguish 
between  the  noise  only  situation  and  a  small  <z2  situation. 

To  use  the  modified  cross  power  spectrum  algorithm  at  low  SXRs  (  i.e.. 
SXR  <  70  ),  we  need  to  investigate  the  contour  plot  some  more. 


SNR=1 0  DELAY=0  «1=a2  =  .23 


Figure  35.  Contour  plot  of  the  modified  cross  power  spectrum  (case  1) 
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Figure  40.  Contour  plot  of  the  modified  cross  power  spectrum  (noise  only). 
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DOPPLER  AXIS 


APPENDIX  E.  MODIFIED  CROSS  POWER  SPECTRUM  PROGRAM 


c* 

c 

c 

c 

c 

c 


•***********Vf***********Vr*iV**i'rVf********************************** 


main  program  of  approach  one. 

this  program  share  subroutine  with 
the  main  program  approach  two 


COMPLEX  XC-128: 255),Y(-128:  255), A(0:  2500), B( 0:  2500) 
COMPLEX  BPF1C-128: 19 1) , BPF2( - 128: 191),BPF(0: 100) 
INTEGER  SNR, BINDEL, DELAY 

DATA  X,Y,A,B/5770*(0.  ,0. )/ ,BPF1 , BPF2 , BPF/741*( 0.  ,0.  )/ 
OPEN ( UN IT=7,FILE=’ FILENAME  FILETYPE  A') 

PI=AC0S( -1. ) 


c - 

c 

* 

c 

input  part 

c 

fs 

-  sampling  frequency 

* 

c 

f 

-  carrier  frequency 

it 

c 

alphal 

-  doppler  effect  at  sensorl 

* 

c 

alpha2 

-  doppler  effect  at  sensor2 

* 

c 

delay 

-  signal  receiving  time  difference 

it 

c 

between  two  sensors 

c 

SNR 

-  signal  to  noise  ratio 

it 

c 

iseed 

-  noise  generator  seed(odd  number) 

it 

c 

c 

n 

-  number  of  data  during  one  second 

* 

FS=64. 

F=23. 

ALPHA 1=. 

23 

ALPHA2=- 

.  02 

DELAY=0 

SNR=100 

N=64 

T=2. *PI* 

F/FS 

ISEED=13 

-At 

c 

c 

input 

data  sampling  at  two  sensors 

it 

it 

c 

* 

c - 

CALL  SI GNAL( X , F+ALPHA 1 , FS , 0 , SNR , I SEED , 0 , 0 ) 

CALL  SIGNAL( Y , F+ALPHA 2 , FS , DELAY , SNR , ISEED ,0,0) 
DO  100  K=-128, 191 
DO  50  J=K,K+63 
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c - * 

c  * 

c  linear  transformation  using  band  pass  filter  * 

c  * 

c - * 

OMEGA=T*FLOAT ( J -K) 

BPF1 (K)=BPF1(K)+X( J)*CMPLX(COS(OMEGA) , -SIN(OMEGA) ) 

50  BPF2(K)=BPF2(K)+Y( J)*CMPLX( COS (OMEGA) , SIN( OMEGA) ) 

c - * 

c  * 

c  normalization  of  BPF2  * 

c  * 

c - * 

100  BPF2(K)=BPF2(K) /CABS( BPF2(K) ) 

DO  300  BINDEL=-128, 127 
DO  200  K=0 , N 

c - * 

c  * 

c  input  of  AR  model  * 

c  * 


200  BPF(K)=BPF1(K)*BPF2(K+BINDEL) 

CALL  CLEAR(A,B,2500) 

CALL  BURGAR( BPF,N+1 , A, IP , VAR) 

CALL  CHANGEFFT(B , A , 11 ,MAX) 

DO  250  I=-32 , - 1 

250  WRITE( 7 , 997 )BINDEL,FLOAT( I )/32. ,CABS(1. /BC 2048+1 ) )*VART 

*,CABS( 1.  / B C 2048+1) ) 

DO  300  1=0,32 

300  WRITE( 7 ,997)BINDEL,FL0AT( I)/32. ,CABS( 1. /B( I) )*VART 

*,CABS( 1. /B ( I) ) 

CL0SE( 7 ) 

STOP 

997  F0RMAT( 1X,I4,2X,E13. 6,2X,E13. 6,2X,E13. 6) 

END 
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APPENDIX  F.  COHERENCE  PROGRAM 


C:V**?V*****yrfr*?VV«V**?V**?VsV5V**#r*******A*5V**********************: 


C 

C 


C 

c*- 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


main  program  of  approach  two 


•>*- «*«  WU  »'#  «!■  • 


■  ^  ^  fct 


COMPLEX  X( - 128:  255), Y( -128: 255), A(0: 2500), B( 0: 2500) 
COMPLEX  BPFlC-128: 191) ,BPF2( -128: 191),BPF(0: 100) 
COMPLEX  AUT0X(0:  100) , AUT0Y( 0:  100),SUMXY 
INTEGER  SNR, BINDEL, DELAY 

DATA  X,Y,A,B/5770*(0.  ,0.  )/ ,BPF1 ,BPF2 , BPF/74l*( 0.  ,0. )/ 
DATA  AUTOX, AUTOY/202*(0. ,0. )/ 

OPEN(UNIT=7,FILE=' FILENAME  FILETYPE  A’) 

PI=AC0S( -1.  ) 


input  part 

fs  -  sampling  frequency 

f  -  carrier  frequency 

alphal  -  doppler  effect  at  sensorl 

alpha2  -  doppler  effect  at  sensor2 

delay  -  signal  receiving  time  difference 
between  two  sensors 
SNR  -  signal  to  noise  ratio 
iseed  -  noise  generator  seed(odd  number) 
n  -  number  of  data  during  one  second 


ALPHA1=. 23 
ALPHA2=-. 02 
SNR=20 


* 


F=23. 

N=64 

FS=64. 

T=2. *PI*F/FS 

ISEED=23 

DELAY=0 


c 

c  input  data  sampling  at  two  sensors  * 

c  * 

c - * 

CALL  S I GNAL( X , F+ALPHA 1 , FS , 0 , SNR , ISEED, 0,0) 

CALL  SIGNAL( Y , F+ALPHA2 , FS , DELAY , SNR , I SEED ,1,0) 

DO  100  K=-128 , 19 1 


DO  50  J=K,K+63 

c - * 

c  * 

c  linear  transformation  using  band  pass  filter  * 

c  * 
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OMEGA=T*FLOAT( J-K) 

BPF1(K)=BPF1(K)+X(J)*CMPLX(C0S( OMEGA), -SIN( OMEGA)) 

50  BPF2(K)=BPF2(K)+Y( J)*CMPLX( COS( OMEGA) ,  SIN( OMEGA) ) 

100  CONTINUE 
SUMX=0. 

DO  150  K=0,N 

c - * 

c  * 

c  fine  AR  order  at  sensorl  * 

c  * 

c - * 

150  AUT0X(K)=BPF1(K) 


CALL  CLEAR(A,B,2500) 

CALL  BURGAR( AUTOX , N+ 1 , A , I P , VARX ) 
CALL  CHANGEFFT( B , A , 1 1 , MAX) 
XMAX=CABS( 1. /B(MAX) ) 

DO  300  BINDEL=-128 , 127 


DO  200  K=0,N 

c - * 

c  * 

c  input  of  AR  model  * 

c  * 

c - * 

AUTO Y ( K ) =BPF2 ( K+B INDEL) 

200  BPF(K)=BPF1(K)*(BPF2(K+BINDEL) ) 

c - - - - - - - - - - -  -  - - * 

c  * 

c  find  auto  power  spectrum  of  sensor2  at  delay  k  * 

c  * 


CALL  CLEAR(A,B, 2500) 

CALL  BURG( AUTOY , N+ 1 , A , I P , VARY ) 

CALL  CHANGEFFT( B , A, 11 ,MAX) 

YMAX=1.  /CABS( B(MAX) ) 

c - - - . - . - . - . * 

c  * 

c  find  auto  power  spectrum  of  cross  term  * 

c  * 

c - * 

CALL  CLEAR( A,B , 2500) 

CALL  BURGCBPF.N+l ,A, IP,VARXY) 

CALL  CHANGEFFT ( B , A , 1 1 , MAX ) 
XSMAXKXMAX*YMAX)**2*VARX*VARY/VARXY*fs 
DO  250  I=-32 , - 1 

250  WRITE(7,997)BINDEL,FLOAT(I)/32. ,CABS(1. /B( 2048+1) )/SQRT(XSMAX)*8. 
* ,CABS( 1.  / B ( 2048+1) ) 

DO  300  1=0,32 

300  WRITE( 7 , 997)BINDEL,FLOAT( I ) / 32. ,CABS(1. /B( I) )/SQRT(XSMAX)*8. 

*,CABS( 1. /B( I) ) 

CL0SE( 7 ) 

STOP 


997  F0RMAT( IX, 14 , 2X,E13.  6 , 2X,E13. 6 , 2X,E13. 6) 
END 
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oooooooooooo  *-■  ooooooonoon 


■>c 


SIGNAL  GENERATER 

INPUT 

N 

-  NUMBER  OF  DATA  POINTS 

FI  ,F2 

-  FREQUENCY  OF  1ST, 2ND  SIGNAL 

AMP1 , 

AMP2 -FREQUENCY  OF  1ST, 2ND  SIGNAL 

OUTPUT 

A(N) 

-  SIGANL 

SUBROUTINE  S I GNAL( B , F , FS , D , SNR , I SEED , ID , NOI SE ) 
COMPLEX  B( -128: 255) 

REAL  A(0: 1) 

INTEGER  SNR,D 
PI=ACOS( - 1.  ) 

DATA  A/0.  ,0.  / 

SIGMA=1. 

A( 0)=SQRT( 2.  *10.  **( SNR/ 1 0 ) ) 

IF(N0ISE. EQ. 1)  A(0)=0. 

DO  100  I=- 128 , 255 
M=( I -D) / 128 

IF(  I.  LT.  D.  OR. M. NE. 0)  M=1 
T=AM0D(F*FL0AT( I-D) ,FS) 

CALL  GAUSS(ISEED, SIGMA ,0. , RANDOM ) 

IF(ID. EQ. 0)  THEN 

X=COS( 2. *PI*T/FS ) *A( M ) +RAND0M 

ELSE 

X=SIN( 2. * p i *t/ F S ) * A ( M ) +RANDOM 
ENDIF 

B( I )=CMPLX(X, 0.  ) 

RETURN 

END 


BURG  ALGORITHM 

* 

INPUT 

* 

N 

NUMBER  OF  DATA  POINTS 

* 

X 

INPUT  SIGNAL 

Vr 

OUTPUT 

* 

IP 

ORDER  OF  AR 

■>v 

A( 0: IP)- 

AR  COEFFICIENTS 

* 

VAR 

DRIVING  NOISE  VARIATION 

Vr 

Vc 


SUBROUTINE  BURGAR(X,N, A, IP, VAR) 

COMPLEX  X(N) , A( 0: N) ,EFK(500) ,EBK(500) 

COMPLEX  EFK1(500) ,EEK1(500) ,AA(20,20) ,SUMN,SUMD 
REAL  RH0( 0: 80) ,FPE( 0: 80) 

INTEGER  START 
RH0( 0)=0. 

FPE( 0)=1.  E64 
A( 0)=CMPLX( 1.  ,0-  ) 

IP=1 


72 


START=1 
DO  10  1=1, N 

10  RHO ( 0 ) =RHO( 0 ) +C AB  S ( X ( I ) ) **2 /FLOAT ( N ) 

DO  20  1=2, N 
EFK1( I)=X( I) 

20  EBK1( I - 1)=X( I - 1) 

LOOP 

K=IP 

SUMN=CMPLX( 0.  ,0.  ) 

SUMD=CMPLX(  0.  ,0.  ) 

DO  30  I=K+1,N 

SUMN=SUMN+EFK 1 ( I )*CONJG( EBKl(I-l)) 

3 0  SUMD=SUMD+C AB S ( EFK 1 ( I ) **2 ) +C AB S ( EBK 1 ( I  - 1 ) **2 ) 

CO  SUMD=SUMD+CABS(EFK1( I ) )**2+CABS( EBK1( I - 1 ) )**2 

AA(K,K)  =  -2.  *SUMN/SUMD 
TEMP=1. -CABS(AA(K,K) )**2 

IF(TEMP. LE. 0. )  TEMP=1. E- 10 
RHO ( K ) =TEMP*RHO( K - 1 ) 

IF(K.  GT.  1)  THEN 
DO  40  J=1 ,K- 1 

40  AA(J,K)=AA(J,K-1)+AA(K,K)*C0NJG(AA(K-J,K-1)) 

END  IF 

DO  60  I=K+2 ,N 

EFK( I)=EFK1( I)+AA(K,K)*£BK1( 1-1) 

60  EBK( I- 1 )=EBK1( I-2)+C0NJG( AA(K,K) )*EFK1( I - 1) 

DO  70  I=K+2 , N 
EFK1( I)=EFK( I) 

70  EBK1(I-1)=EBK(I-1) 

IF(N-K. EQ. 1)  THEN 

FPE(K)=FPE(K-1)+1. 

ELSE 

FPL  ( K )  =RHO(  K )  '-''FLOAT ( N+ 1+K )  /  FLOAT ( N  - 1  -  K  ) 

END  IF 
IP=IP+1 

UNTIL(  FPE(K).  GT.  FPE(K-l).  AND.  K.  GT. START) 

IP=K- 1 

DO  100  1=1, IP 
100  A( I )=AA( I , IP) 

VAR=RHO( IP) 

RETURN 

END 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BURG  ALGORITHM 

Vr 

INPUT 

* 

N  -  NUMBER  OF  DATA  POINTS 

* 

X  -  INPUT  SIGNAL 

* 

IP  -  ORDER  OF  AR 

?v 

OUTPUT 

* 

A( 0:  IP) -  AR  COEFFICIENTS 

Vf 

VAR  -  DRIVING  NOISE  VARIATION 

‘V 

* 
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SUBROUTINE  BURG(X,N,A, IP, VAR) 

COMPLEX  X(N) , A( 0:  N) ,EFK(500) ,EBK(500) 

COMPLEX  EFK1(500) ,EBK1(500) ,AA(20,20) ,SUMN,SUMD 
REAL  RH0(0: 80) ,FPE(0: 80) 

INTEGER  START 
RHO( 0)=0. 

A(0)=CMPLX( 1.  ,0. ) 

DO  10  1=1, N 

10  RHO( 0 ) =RHO( 0 ) +C ABS ( X( I ) ) **2 / FLOAT( N ) 

DO  20  1=2, N 
EFK1( I)=X( I) 

20  EBK1( I-1)=X( 1-1) 

DO  1000  K=1 , IP 

SUMN=CMPLX( 0. ,0. ) 

SUMD=CMPLX( 0. ,0. ) 

DO  30  I=K+1,N 

SUMN=SUMN+EFK 1 ( I ) *CON JG ( E  BK 1 ( I - 1 ) ) 

3  0  SUMD=SUMD+CAB S ( EFK 1 ( I ) ) **2+C  AB  S ( EBK 1 ( I - 1 ) ) **  2 

AA(K,K)=-2.  *SUMN/SUMD 
TEMP=1.  -CABS( AA(K,K) )**2 

IF(TEMP.  LE.  0.  )  TEMP=1.  E-10 
RHO ( K ) =TEMP*RHO ( K - 1 ) 

IF(K.GT.  1)  THEN 
DO  40  J=1,K-1 

40  AA( J,K)=AA( J,K-1)+AA(K,K)*C0NJG(AA(K-J,K-1) ) 

END  IF 

DO  60  I=K+2,N 

EFK( I )=EFK1( I )+AA(K ,K)*EBK1( I -1) 

60  EBK( I -1)=EBK1( I-2)+C0NJG( AA(K,K) )*EFK1( I -1) 

DO  70  I=K+2 ,N 
EFK1( I)=EFK( I) 

70  EBK1(I-1)=EBK(I-1) 

1000  CONTINUE 

DO  100  1=1, IP 
100  A( I )=AA( I , IP) 

VAR=RHO( IP) 

RETURN 

END 


* 

CHANGE  FFT 

* 

INPUT 

* 

B 

-  AR  COEFFICIENTS 

Vr 

A 

-  TEMPORARY  USING  ARRAY 

Vr 

ISIZE 

-  ALOG(FFT  DATA  POINTS) /AL0G2 

Vr 

OUTPUT 

Vf 

A 

-  POWER  SPECTRUM 

Vr 

Vr 
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SUBROUTINE  CHANGEFFT( A ,  B ,  I S I ZE ,  K  ) 
COMPLEX  A( 0 : 2**ISIZE-1) ,B(0: 2**ISIZE-1) 
INTEGER  REVRSE 
PI=AC0S( -1.  ) 

N=2**ISI2E-1 

K=0 

FS=FL0AT(N+1) 

CALL  FFT(N, ISIZE , A , B ,FS) 

CALL  FINDMAX(A,N,K) 

RETURN 

END 


Vr 


FINDMAX 

* 

INPUT 

* 

N 

-  NUMBER  OF  SPECTRUM  POINTS(N+l) 

Vr 

A 

-  POWER  SPECTRUM 

OUTPUT 

* 

MAX 

-  ARRAY  INDEX  OF  MAXIMUM  POWER  SPECTRUM 

■jV 

oo 


SUBROUTINE  FINDMAX( A , N , MAX) 
COMPLEX  A( 0: N) 

MAX=0 


AMAX=1. E66 
DO  100  1=0, N 

IF(CABS(A( I) ). LT. AMAX)  THEN 
MAX=I 


AMAX=CABS( A( I ) ) 

END  IF 
CONTINUE 
RETURN 
END 


REVERSE  * 

BIT  REVERSE  ORDER  INDEX  CHANGING  SUBROUTINE 
INPUT  * 

N  -  ARRAY  INDEX  * 

ISIZE  -  ALOG(FFT  DATA  POINTS)/ALOG2 

OUTPUT  * 

REVERSE-  BIT  REVERSE  ORDER  INDEX  * 

* 


INTEGER  FUNCTION  REVRSE(N, ISIZE) 
INTEGER  A( 20) 

DO  200  1=1, ISIZE 
A(I)=MOD(N,2) 

N=N/2 
200  CONTINUE 

REVRSE=0 
DO  300  1=1 , ISIZE 
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REVRSE=REVRSE+A( I)*2**( ISIZE-I ) 
300  CONTINUE 
RETURN 
END 


C . 

c 

C  DFT 

C  COOLEY  -  TUKEY  ALGORITHM 


c 

INPUT 

* 

c 

N 

-  NUMBER  OF  DATA  POINTS 

-V 

c 

ISIZE 

-  ALOG(N)/ALOG2 

■sV 

c 

A 

-  BIT  REVERSE  ORDER  INDEXED  TIME  DOMAIN 

Vc 

c 

B 

-  TEMPORARY  ARRAY 

Vc 

c 

OUTPUT 

* 

c 

A 

-  FREQUENCY  DOMAIN 

»v 

c 

Vc 

C . 

-  -Vf 

SUBROUTINE  DFT(N, ISIZE , A,B) 

COMPLEX  A( 0 :  N) , B( 0: N) ,W 
PI=AC0S( -1. ) 

DO  500  11=1, ISIZE 
I SAGE 1=2**( I 1 - 1 ) 

ISTAGE=2*'T1 
DO  200  1=0, N 

ITEST=MOD( I , ISTAGE) 

IF( ITEST.  GT.  ISAGE1)  THEN 
K=ITEST-ISAGE1 

THEATA=2. *PI*FLOAT(K)/FLOAT( ISTAGE) 
T1=SIN(THEATA) 

T2=COS(THEATA) 

IF(ABS(T1). LT.  IE-5)  T1=0 
IF(ABS(T2). LT.  IE-5)  T2=0 
W=CMPLX(T2,-T1) 

A( I )=A( I )*W 

END  IF 

200  CONTINUE 

DO  300  1=0, N 

ITEST=MOD( I , ISTAGE) 

IF( ITEST. GE. ISAGE1)  THEN 

B( I)=-A( I)+A( I -ISAGE1) 

ELSE 

B( I)=A( I)+A( I+ISAGE1) 

END  IF 

300  CONTINUE 

DO  400  1=0, N 
A( I )=B( I) 

400  CONTINUE 

500  CONTINUE 
RETURN 
END 
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C-... 

c 

' 

c 

FFT 

* 

c 

1. INDEX  CHANGE 

* 

c 

2. USE  COOLEY  TUKEY  ALGORITHM 

iV 

c 

INPUT 

it 

c 

N 

-  NUMBER  OF  DATA  POINTS 

* 

c 

ISIZE 

-  AL0G(N)/AL0G2 

* 

c 

B 

-  INPUTCTIME  DOMAIN) 

* 

c 

FS 

-  SAMPLING  FREQUENCY 

* 

c 

OUTPUT 

if 

c 

A 

-  FREQUENCY  DOMAIN 

if 

c 

if 

p.  .  . . 

SUBROUTINE  FFT(N, ISIZE 
COMPLEX  A(0: N) ,B(0:  N) 
INTEGER  REVRSE 

DO  100  1=0, N 

,A,B ,FS) 

K=I 

100 

A( REVRSE (K,ISIZE))= 

B(  I) 

CALL  DFT(N, ISIZE , A, 
DO  200  1=0, N 

B) 

200 

A( I )=A( I )/FS 

RETURN 

END 

p _ 

_  _ _ ...» _ * 

c 

it 

c 

CLEAR 

* 

c 

INPUT 

5f 

c 

N 

-  NUMBER  OF  DATA  POINTS 

* 

c 

A,B 

-  INPUT  arrays 

it 

c 

OUTPUT 

if 

c 

A ,  B 

-  reinitiallized  arrays 

if 

c 

if 

C---. 

SUBROUTINE  CLEAR(A,B,N) 

COMPLEX  A(0:N),B(0:N) 

DO  100  1=0, N 

100  A( I)=B( I)=CMPLX( 0.  ,0.  ) 

RETURN 
END 
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